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ABSTRACT 

First principles linear response calculations are used to investigate the lattice dynamics of 
what is thought to be the third most abundant phase in the lower mantle, CaSiOa 
perovskite. The commonly assumed cubic structure (Pm3m) is found to be dynamically 
unstable at all pressures, exhibiting unstable modes along the Brillouin zone edge from the 
M-point to the R-point. Based on these results, we predict that the ground state structure 
of CaSi0 3 perovskite is a distorted phase with lower than cubic symmetry. We predict 
that a phase transition occurs in CaSiOs perovskite within the earth's lower mantle from 
the low temperature distorted phase to the cubic phase at high temperature. The 
predicted phase transition provides a possible explanation of some of the seismological 
observations of reflective features within the lower mantle. 
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INTRODUCTION 



CaSi0 3 perovskite is thought to comprise between 6 and 12 weight % of the lower half 
of the earth's transition zone and lower mantle (depths between 500 and 2900 km) (Irifune 
1994; O'Neill and Jeanloz 1990; Ita and Stixrude 1992). Its structure throughout this 
regime is generally assumed to be cubic because x-ray diffraction studies have found no 
detectable deviation from Pm3m symmetry (Liu and Ringwood 1975; Mao et al. 1989; 
Wang et al. 1996). Theoretical studies based on ionic models (Wolf and Jeanloz 1985; 
Hemley et al. 1987; Wolf and Bukowinski 1992), periodic Hartree-Fock (Sherman 1993), 
and pseudopotential calculations (Wentzcovitch et al. 1995) have supported this picture. 

Here, we go beyond previous studies by investigating from first principles the full 
phonon spectrum of cubic CaSi0 3 perovskite from low pressures to those typical of the 
lower mantle. The calculations are based on the Linearized Augmented Plane Wave 
(LAPW) method, widely accepted as among the most accurate methods for solving the 
band structure/total energy problem. We find that the cubic phase is dynamically 
unstable at all pressures and that the ground state structure must possess lower than cubic 
symmetry. 

These results are important for our understanding of the physics and chemistry of this 
mantle phase - its phase stability, elasticity, and its ability to incorporate other cations 
such as Mg, Fe, and Al, for example, will be affected by its symmetry. Moreover, the 
existence of a low symmetry ground state raises the possibility of a temperature induced 
phase transition in this mineral within the earth. Phase transitions in lower mantle 
constituents have important implications for our understanding of recent seismological 
observations of reflective features within the lower mantle (Revenaugh and Jordan 1991; 
Kawakatsu and Niu 1994; LeStunff et al. 1995) near 700, 900, and 1200 km depth. These 
observations indicate the presence of previously unpredicted phase transitions or sharply 
bounded compositional heterogeneities, challenging our traditional view of this region as 
being compositionally and mineralogically homogeneous. 
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COMPUTATIONS 



The LAPW linear response calculations determine the self-consistent first order 
response of the electronic charge density to generalized perturbations (Yu and Krakauer 
1994). The calculations are based on density functional theory; the only essential 
approximation is to the exchange-correlation potential. We use the well studied local 
density approximation (LDA), which has been applied successfully in studies of silicates, 
including Mg-perovskite (Stixrude and Cohen 1993). We compute the linear response to 
shifts in nuclear positions and imposed external fields, which yield the elements of the 
dynamical matrix, the dielectric constant (e) and the Born effective charges (Z*), from 
which the phonon frequencies are determined (Yu and Krakauer 1995; Lee and Gonze 
1994). The perturbations need not be commensurate with the unit cell, so that points in 
the Brillouin zone away from the zone center are readily investigated. We calculate the 
dynamical matrix of the cubic Pm3m phase of CaSi0 3 perovskite at four points in the 
Brillouin zone (r, X, M, R). Computational variables (/c-point mesh: 4x4x4 for 
computations of the dynamical matrix elements, 8x8x8 for e and Z*; number of basis 
functions per atom 150) were chosen such that phonon frequencies are converged to 
better than one wavenumber. The full phonon dispersion curve is determined using an 
interpolation scheme that separates short-range forces from long-range, Coulombic 
interactions (Yu and Krakauer 1994; Gonze et al. 1994). 

RESULTS 

We find that the Born effective charges in CaSiOs perovskite differ significantly from 
formal ionic charges; by as much as le (Fig. 1). Compression causes Z* of Si and Ca to 
deviate farther from their formal ionic values, and Z* of O to approach -2. The dielectric 
constant, e = 3.8 at zero pressure, lies between the known values for CaO (Lide and 
Frederikse 1994) and Si0 2 stishovite (Stishov and Popova 1961) and varies slowly with 
pressure (Fig. 1). LO-TO splitting of zone-center vibrational frequencies is on the order of 
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50-200 wavenumbers (Fig. 2). All zone-center modes are found to be dynamically stable. 
An examination of the eigenvectors shows that the highest frequency modes involve 
stretching of the octahedral Si-0 bond. The change in frequency of these modes with 
compression is greatest because the length of the Si-0 bond shrinks in direct proportion to 
compression in the cubic structure. 

The full phonon dispersion curves reveal dynamical instabilities in the cubic structure 
along the zone boundaries (Fig. 2). Instabilities occur at the M- and R-points and along 
the zone edges from M to R. The corresponding imaginary frequencies grow in magnitude 
upon compression. An examination of the eigenvectors associated with these unstable 
modes shows that they consist of coupled rotations of the Si0 6 octahedra. The unstable 
modes at the M-point (M 2 ) and R-point (R25) are associated, respectively, with in-phase 
and out-of-phase rotations of the octahedra about [100]. Instabilities along M-R were 
previously found in ionic model calculations, but only at pressures (P) above 80 GPa 
(Hemley et al. 1987; Wolf and Bukowinski 1992). Here we find that the instabilities are 
more profound and that they persist even in the expanded lattice (P—-8 GPa). Unstable 
modes of the same type have also been found in theoretical investigations of cubic MgSi0 3 
perovskite which account for the observed orthorhombic (Pbnm) symmetry of this 
material (Wolf and Jeanloz 1985; Hemley et al. 1987; Wolf and Bukowinski 1992). 

Our results indicate that the cubic Pm3m structure of CaSi0 3 perovskite is 
dynamically unstable and that the ground state of this material must possess lower 
symmetry. However, the ground state is not of direct relevance to the lower mantle where 
temperatures exceed 2000 K. Even if the cubic phase is dynamically unstable, it may be 
thermodynamically stable at high temperature because of its greater entropy (e.g. Salje 
1990). 

The temperature at which the low temperature distorted phase transforms to the high 
temperature cubic phase depends not only on the unstable mode frequencies, but also on 
the energetics of finite displacements along the unstable mode eigenvectors. Using LAPW 
total energy calculations, we find the total energy as a function of displacement along the 
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most unstable mode eigenvector (R25) and the associated minimum energy displacement 
[frozen phonon approach; Cohen 1992]. The symmetry of the structure associated with this 
finite displacement is tetragonal I4/mcm, with 10 atoms in the unit cell. The results at 
mid-mantle pressures (P=80 GPa) show that the minimum energy displacement 
corresponds to an octahedral rotation angle of 7 degrees and is 360 K per octahedron lower 
in energy than the cubic phase (Fig. 3). The relatively small minimum energy octahedral 
rotation angle may explain why previous x-ray diffraction experiments have failed to detect 
deviations from cubic symmetry. Assuming that the octahedral rotations are rigid 
(O'Keefe et al. 1979), the corresponding c/a ratio of the distorted phase differs by only 0.7 
% from the cubic structure - below the detection limit of previous non-hydrostatic (Mao et 
al. 1989) and quasihydrostatic (Wang et al. 1996) experiments (H. K. Mao, personal 
communication) . 

These results constrain the parameters of a simple model which we use to estimate the 
transition temperature from the distorted 1 4 /mem structure to the cubic Pm3m structure. 
The model (Bruce 1980) consists of on-site terms which correspond to finite displacements 
along the unstable zone-boundary mode eigenvectors, and inter-site couplings 

T ran 

U = EV(Q t ) + iE(Q 3 + Q l ) 2 (i) 

where U is the energy contained in octahedral rotations, Q is the normal coordinate of the 
octahedral rotation, and V(Q) is the on-site term given by our total energy results 
(Fig. 3). The form of the nearest-neighbor inter-site coupling term is determined by the 
requirement that it vanish for a pure R25 mode distortion of the crystal, in which 
neighboring octahedra rotate in opposite directions. The intersite coupling parameter, J, 
is most simply given by the curvature along T-R of the unstable mode eigenenergy 
evaluated at the R-point (Fig. 2). The relative magnitudes of the on-site and inter-site 
coupling terms are such that the phase transition is expected to occur in the displacive 
limit. In this case, if we ignore coupling between the unstable mode eigenvector and other 
modes, the transition temperature is given by 
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where k B is the Boltzmann constant, and g(3) = 0.5054. The estimated transition 
temperature is similar to estimates of temperatures in the lower mantle at similar 
pressures: P=80 GPa corresponds to a depth of 1850 km and a temperature of 
approximately 2500-3000 K. 

DISCUSSION 

Our estimate of T c may be uncertain by several hundred K and is likely to be a lower 
bound since we have ignored coupling between unstable mode eigenvectors - this may lower 
the total energy of the distorted phase further relative to the cubic phase. Nevertheless, 
the fact that our estimated T c is comparable to lower mantle temperatures indicates that a 
phase transition in CaSiOs perovskite is likely to occur in the lower mantle. 

The predicted phase transition in CaSi0 3 perovskite is expected to be associated with 
an elastic anomaly which may be the cause of at least a subset of the reflective features 
near 700, 900, and 1200 km depth in the lower mantle (Revenaugh and Jordan 1991; 
Kawakatsu and Niu 1994; LeStunff 1995). Other phase transitions have been discussed as 
possibly occurring in the earth's lower mantle. A phase transition in the lower mantle's 
most abundant constituent, MgSi0 3 perovskite, has been suggested (Meade et al. 1995), 
but theoretical computations show such a transition unlikely for reasonable geotherms 
(Stixrude and Cohen 1993; Warren and Ackland 1996). Recent results show that Si0 2 
undergoes a phase transition from the stishovite to the CaCl 2 structure at lower mantle 
pressures (Cohen 1992; Kingma et al. 1995). Our results on CaSiOs perovskite indicate 
the presence of yet another phase transition within the lower mantle, the first high 
pressure phase transition predicted with the linear response method. 
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Fig. 1. Born effective charges (left axis) and dielectric constant (right axis) as a 
function of compression. Pressures range from -8 GPa (V—310 Bohr 3 ) to 140 GPa (V=220 
Bohr 3 ) and were determined by fitting a third order finite strain expression to the LAPW 
total energy as a function of volume. The two symmetrically distinct elements of the 
oxygen Born effective charge tensor are shown, Z Ql for motion of O along the Si-0 bond, 
and Zo 2 for motion normal to the bond. Dashed lines indicate formal ionic charges. 

Fig. 2. Phonon spectrum of the cubic phase at V=310 Bohr 3 (P=-8 GPa) (top) 
compared with that at V=2A0 Bohr 3 (P=80 GPa) (bottom). All zone-center modes 
(T-point) are stable. However, zone-boundary modes at the M- and R-points are unstable, 
as shown by the existence of imaginary frequencies. The final panel (M-R) shows that the 
edges of the cubic Brillouin zone are everywhere unstable. The magnitude of the imaginary 
eigenfrequencies increases with compression. Symmetry designations are those of Cowley 
(1964) as corrected by Boyer and Hardy (1981). At T there are four TO modes (r 5 and 
inactive T 2 5 symmetries) and three LO modes (Pi). The volume dependence of the zone 
center modes is described by u = co (V /Vq)"' ', where the frequencies at V = 310 Bohr 3 (near 
zero pressure), u , are: 197, 274, 327, 393, 621, 691, 926 cm -1 ; and the respective mode 
Gruneisen parameters, 7$, are: 2.5, 1.5, 1.2, 1.2, 0.8, 2.0, 1.2. 



11 



Fig. 3. Total energy (per octahedron) of the I A/ mem tetragonal structure as a 
function of finite displacement along the R 25 eigenvector at V=240 Bohr 3 (P=80 GPa). 
The LAPW total energy results (symbols) are fit to a polynomial in the square of the 
normal coordinate, Q: V(Q) = AQ 2 /2 + BQ 4 /A. The normal coordinate per octahedron, 

Q = I m i x h where the sum is over the cartesian displacements (x i: in Bohr) of the n 

1=1 

atoms in the unit cell, with masses m^, in amu. This expression, with A=-2.04 mRy amu -1 
Bohr -2 , and 5=0.456 mRy amu -1 Bohr -4 constrains the on-site term in the simple model 
of the phase transition discussed in the text. The remaining parameter, J— 1.19 mRy 
amu -1 Bohr -2 , is determined from the dispersion of the unstable mode eigenvalue (Fig. 2). 
The unstable mode frequency in the cubic phase (Q = 0) calculated from the curvature of 
the polynomial fit (173i cm -1 , % = \^T) agrees with the linear response result (164i cm -1 ). 
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